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Abstract 

Radial  basis  function  interpolation  has  an  advantage  over  other  methods  in  that  the 
interpolation  matrix  is  nonsingular  under  very  weak  conditions  on  the  location  of  the 
interpolation  points.  However,  we  show  that  point,  location  can  have  a significant  effect 
on  the  performance  of  an  approximation  in  certain  cases.  Specifically,  we  consider  multi- 
quadric and  thin  plate  spline  interpolation  to  small  data  sets  where  derivative  estimates 
are  required.  Approximations  of  this  type  are  important  in  the  motion  of  unsteady  in- 
terfaces in  fluid  dynamics.  For  data  points  in  the  plane,  it  is  shown  that  interpolation  to 
data  on  a circle  can  be  related  to  the  polynomial  case.  For  scattered  data  on  the  sphere, 
a comparison  is  made  with  the  results  of  Sloan  and  Womersley. 


1 Introduction 

Radial  basis  functions  (RBFs)  such  as  multiquadrics  or  thin  plate  splines  have  been 
successfuly  used  for  scattered  data  approximation  in  many  applications.  They  have  been 
shown  to  perform  well  for  data  fitting,  although  problems  of  ill-conditioning  and  the 
computational  cost  of  processing  large  data  sets  must  be  handled  carefully.  In  general, 
when  considering  the  accuracy  of  a RBF  interpolant,  a balance  must  be  achieved  between 
the  reduction  in  fill  distance  necessary  for  convergence  of  the  approximation  to  an  as- 
sumed underlying  function  and  the  need  to  maximise  the  separation  distance  between 
data  points  to  avoid  problems  of  ill-conditioning  [4]. 

In  the  present  study,  we  focus  on  the  use  of  RBF  approximation  as  one  stage  of  a 
larger  algorithm  to  compute  the  evolution  of  an  unsteady  interface  in  fluid  dynamics. 
The  accuracy  of  the  approximations  made  in  the  algorithm  and  the  interaction  between 
its  different  stages  determine  whether  the  output,  is  close  to  the  true  solution  of  the 
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governing  equations  or  whether  spurious  effects  are  produced.  In  the  three-dimensional 
setting,  a typical  example  is  described  by  Zinchenko  et  al.  [8]  where  the  deformation  of 
liquid  drops  in  a viscous  medium  is  studied.  A critical  feature  of  the  algorithm  is  the 
approximation  of  the  normal  directions  and  curvatures  of  the  droplet  surface  defined  at 
a number  of  discrete  points. 

The  focus  here  is  algorithmic  rather  than  theoretical  and  we  investigate  the  perform- 
ance of  multiquadric  and  thin  plate  spline  local  interpolants  applied  to  the  determination 
of  normal  directions  and  curvatures  of  a smooth,  closed  surface.  Certain  configurations 
of  data  points,  such  as  points  located  on  a circle,  impose  constraints  on  the  interpolant. 
A framework  for  understanding  the  behaviour  of  the  RBF  interpolants  is  provided  by  a 
comparison  with  the  multivariate  polynomial  interpolant  of  de  Boor  and  Ron  [1]  and  by 
considering  the  free  parameter  in  the  multiquadric  as  a tensioning  parameter  [2]. 

2 Approximation  method 

A common  approach  to  solving  fluid  dynamics  problems  that  include  moving  inter- 
faces, combines  a computational  grid  with  meshless  approximation  methods.  The  gov- 
erning partial  differential  equations,  or  corresponding  integral  equation  formulation,  are 
solved  on  the  grid,  while  quantities  characterising  the  interface  are  computed  as  meshless 
scattered  data  approximations. 

Here  we  examine  the  behaviour  of  local  RBF  approximations  in  the  general  context 
described  by  Zinchenko  et  al.  [8].  For  a given  data  set,  a particular  point  is  selected 
together  with  its  nearest  neighbours  giving  a set  of  typically  6 or  7 points.  The  initial 
locations  of  these  points  may  be  determined  by  a regular  mesh,  but  the  surface  is  allowed 
to  deform  so  that  the  approximation  is  essentially  to  a small  set  of  scattered  data.  The 
constructed  RBF  interpolant,  S,  can  be  expressed  as 

N K 

S(x)  = °i^(  I lX  _ Xi  1 1 ) + Yj  biPi  (*): » 

j— 1 *=1 

with  the  constraint 

N 

ajPi(xj)  = 0,  for  1 < i < K, 

3= i 

where  x 6 if?2  and  {pi(z)}i=i:/c  is  a basis  for  the  space  of  bivariate  polynomials  of  degree 
< m — 1 with  K = m(m  + 1 ) /2 . The  chosen  forms  for  <j)  are  the  thin  plate  spline 

^(llz-abll)  = ||x  — Xjl^logUz-OJjll,  (TPS) 

and  the  multiquadric 

4>{\\x  ~ Xj\\)  = {\\x-Xj\\2 + <?)$,  (MQ) 

with  1 1 • 1 1 taken  to  be  the  Euclidean  norm. 

A framework  for  interpreting  the  computed  results  in  the  context  being  considered 
can  be  derived  from  [2]  where  the  arbitrary  parameter,  c,  of  the  MQ  function  is  viewed 
as  a tensioning  parameter.  As  c — > oo  the  MQ  interpolant  approaches  the  correspond- 
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ing  polynomial  interpolant  to  the  given  data,  while  as  c — + 0,  the  MQ  surface  is  ten- 
sioned. Multivariate  polynomial  interpolation  can  fail  on  particular  point  sets  and  this 
has  provided  a motivation  for  using  RBF  methods.  However,  the  algorithm  of  de  Boor 
and  Ron  [1]  provides  a reliable  means  of  computing  the  ‘least’  polynomial  interpolant. 
This  algorithm  is  used  to  compute  a polynomial  fit  as  one  reference  point  for  the  in- 
terpretation of  the  MQ  interpolants.  A second  reference  point  is  provided  by  the  TPS 
interpolant  which  gives  a minimum  energy  surface  in  a certain  norm.  This  is  shown  to 
correspond  closely  to  the  MQ  fit  for  a ‘small’,  but  nonzero  value  of  c.  The  MQ  inter- 
polant can  thus  be  shown  to  connect  the  minimum  energy,  tensioned,  TPS  surface  with 
the  polynomial  fit  to  given  data  as  c increases.  In  a fluid  dynamics  context  a fluid-fluid 
interface  is  often  assumed  to  be  represented  by  a C°°  function  (although  cusps  may 
occur  requiring  a change  in  the  representation).  This  would  suggest  that  a high  degree 
polynomial  would  be  preferred  to  a TPS  surface. 


Fig.  1.  Interpolants  to  random  data  at  6 points  (+)  in  the  plane:  (left)  polynomial 
(upper)  and  multiquadric  (c  = 10)  (lower),  contours  [0:0.1:2] ; (right)  thin  plate  spline 
(upper)  and  multiquadric  (c  = 0.4)  (lower),  contours  [0:0.1:1.1]. 


3 Scattered  data  in  the  plane 

To  illustrate  the  behaviour  of  local  interpolation  by  MQ  and  TPS  methods,  random 
points  in  the  xy-plane  (with  - 1 < x, , ?/,  < 1 , for  i = 1 : 6)  arc  associated  with  random 
data  values,  /*  (-1  < fi  < 1).  Figure  1 shows,  in  the  upper  frames,  the  two  reference 
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Fig.  2.  Effect  of  varying  the  parameter  c on  multiquadric  interpolants  to  random  data 
in  the  plane:  (left)  norms  of  the  difference  between  multiquadric  and  polynomial  in- 
terpolants (upper  curve  doo  = ll-lloo,  lower  curve  d2  — IMb/'v/iV);  (right)  curvature 
(k  = 2 H)  computed  at  the  centroid:  — multiquadric;  — • — thin  plate  spline;  — poly- 
nomial. 

interpolating  surfaces:  (left)  the  polynomial  surface  computed  by  the  algorithm  of  [1] 
and  (right)  the  TPS  surface.  The  lower  frames  give  the  contours  of  the  MQ  interpolants 
for  c = 10.0  (left)  and  c = 0.4  (right).  There  is  a close  correspondence  between  the  upper 
and  lower  frames  on  each  side,  but  a large  difference  between  the  polynomial  and  TPS 
surfaces. 

Figure  2 (left)  shows  the  difference  between  the  MQ  surface  and  the  polynomial 
reference  interpolant  computed  on  a regular  grid  on  the  interior  of  the  circle  with  centre 
at  the  centroid  of  the  data  points  (0.44,  —0.09)  and  radius  the  maximum  distance  from 
the  centroid  to  a data  point.  There  is  convergence  of  the  MQ  surface  to  the  polynomial 
as  1/c  — > 0,  but  the  condition  number  of  the  interpolation  matrix  increases  until  the 
calculation  cannot  be  continued.  For  c = 10.0  the  condition  number  is  3 x 107. 

As  an  indication  of  the  behaviour  of  first  and  second  partial  derivatives  of  the  inter- 
polating surfaces  we  calculate  the  curvature  at  the  centroid  of  the  data  points  for  the 
polynomial  and  TPS,  together  with  MQ  as  c varies,  using  k — 2 H where  H is  the  mean 
curvature.  Figure  2 (right)  shows  that  kmq  for  the  MQ  interpolant  coincides  with  the 
value  ktps  = -0.46  for  the  TPS  when  c « 0.4.  When  c < 0.4,  kMq  < ktps , while 
kmq  —■ ► Kp  = -9.78,  the  polynomial  curvature,  as  c increases. 

An  interesting  example  is  presented  in  [1]  of  polynomial  interpolation  for  points 
located  at  the  vertices  of  a regular  hexagon 


with  data  values  /j  = (— 1)*.  This  gives  the  interpolant 

p(x,  y)  = x3  - 3xy2. 

Since  the  points  lie  on  the  unit  circle,  the  quadratic  polynomial 

p2{x,y)  = 1 -x2  -y2 


(3.2) 
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vanishes  at  the  data  points  and  this  causes  difficulties  for  general  polynomial  methods. 
MQ  interpolants  do  not  suffer  from  these  difficulties.  When  c = 10.0,  the  MQ  surface 
is  very  close  to  (3.2).  As  c becomes  smaller,  the  MQ  surface  approaches  that  of  the 
TPS  with  the  data  values  becoming  local  maxima  or  minima  as  the  surface  is  tensioned. 
In  addition,  the  restriction  of  the  data  points  to  a circle  implies  that  the  interpolating 
polynomial  is  harmonic,  but  the  convergence  of  the  approximation  is  only  first  order  [1]. 
The  MQ  surface  for  large  c inherits  the  properties  of  the  polynomial  fit.  Thus,  points  on 
a circle  are  ‘good’  if  the  data  being  interpolated  correspond  to  a harmonic  function,  but 
‘bad’  if  the  data  describe  a function  which  has  a maximum  or  minimum  within  the  circle 
or  a singularity.  These  constraints  on  the  interpolant  are  discussed  further  in  Section  5. 

4 Scattered  data  on  the  sphere 

In  this  section  we  examine  the  accuracy  obtained  from  three  separate  methods  for  in- 
terpolating scattered  data  on  the  unit  sphere  S2  C SR3.  In  particular  we  compare  the 
results  obtained  using  the  MQ  basis  function  in  9?3  with  those  obtained  using  the  spher- 
ical harmonics  of  Sloan  and  Womersley1  [6]  and  the  C1  Hermite  interpolant  of  Renka  [5]. 
For  the  multiquadric  function,  we  list  the  uniform  norm  interpolation  errors  calculated 
using  a range  of  values  for  the  shape  parameter  c. 


os* . 
o»«- 

OM 

OM 

<>»; 
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Fig.  3.  Minimum  energy  points  and  spherical  cap. 

The  point  distribution  used  is  the  256  ‘minimum  energy’  points  of  Fliege  and  Maier 
[3]  and  the  uniform  norm  interpolation  errors  are  calculated  at  points  distributed  on  a 
spherical  cap  (see  [5]). 

The  following  functions  are  used  for  the  comparisons  in  Table  4,  where  the  results 
presented  in  [7]  are  labelled  ‘W&S’. 

FI  = ^ex+v+z,  F2  = — 5sin(l  + lOz), 

F3=  IMIi/lO,  F4  = sin2(l  + 11x110/10. 

We  note  from  Table  4 that  the  multiquadric  function  provides  consistently  better 
interpolants  to  the  four  test  functions  compared  with  the  spherical  harmonics.  Here,  the 

1Uniform  norm  errors  used  for  comparison  are  approximate  only  and  were  taken  from  graphical  representations 

presented  in  Womersley  and  Sloan  [7]. 
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Method 

FI 

F2 

F3 

F4 

W & S 

2.0000e-10 

0.5000 

0.1100 

0.0500 

Renka 

0.0013 

0.1951 

0.0054 

0.0055 

MQ  c = 0.01 

6.0128e-04 

0.3276 

0.0051 

0.0051 

MQ  c = 1 

4.5807e-10 

0.0175 

0.0076 

0.0062 

MQ  c = 2 

2.2615e-13 

0.0227 

0.0079 

0.0065 

Tab.  1.  Comparison  of  uniform  norm  errors. 


points  have  been  chosen  to  minimise  the  interpolation  errors  for  the  harmonic  functions, 
yet  we  see  from  results  given  in  [7]  that  increasing  the  number  of  points  in  the  distribu- 
tion (which  also  increases  the  degree  of  the  interpolating  function)  does  not  necessarily 
produce  better  accuracy.  However,  these  point  distributions  when  used  for  the  multi- 
quadric function  provide  consistently  better  accuracy.  Further  evidence  suggests  that 
points  considered  optimal  for  the  spherical  harmonics  are  also  ‘good’  for  the  multiquad- 
ric function  when  compared  to  an  equal  number  of  generally  scattered  points.  However, 
this  is  due  to  the  uniformity  of  the  point  distributions  and  similar  results  can  be  obtained 
on  a refined  icosahedral  mesh. 


Method 

12  pts 

92  pts 

362  pts 

Renka 

0.1730 

0.0103 

0.8230e-03 

MQ  c = 0.01 

0.2596 

0.0170 

0.0020 

MQc=l 

0.0715 

7.7662e-05 

1.9678e-10 

MQ  c = 2 

0.0442 

3.8206e-05 

3.4113e-ll 

Tab.  2.  Multiquadric  vs  Renka  for  f(x,y,z)  = sin  (a:  + y)  + sin(:r2)  . 

The  Renka  algorithm  produced  similar  results  to  those  obtained  using  the  multiquad- 
ric (for  small  c)  for  the  F3  and  F4  functions,  although  the  results  for  the  functions  FI  and 
F2  were  poor.  Further  comparisons  with  the  Renka  algorithm  have  been  made  using  12, 
92  and  362  icosahedral  points  to  interpolate  the  function  f(x,y,z ) = sin(x-l-y)  +sin(:E2). 
The  uniform  norm  interpolation  errors  have  been  calculated  on  the  previously  mentioned 
spherical  cap.  Again  we  see  that  the  multiquadric  function  produces  better  accuracy  than 
the  Renka  method  when  the  number  of  interpolation  points  is  increased. 

5 Evolution  of  a smooth  closed  surface 

In  this  section  we  return  to  the  local  interpolation  scheme  of  §3  and  apply  it  to  scattered 
data  on  a smooth  closed  surface.  This  is  the  setting  described  in  [8],  where  initially  the 
interface  is  spherical  with  the  point  locations  determined  by  subdivision  of  an  icosahedral 
mesh.  Each  set  of  points  consists  of  a central  point  together  its  nearest  neighbours, 
giving  sets  of  6 points  associated  with  the  12  vertices  of  the  icosahedron  and  sets  of 
7 points  otherwise.  The  local  method  of  Renka  [5]  is  followed  and,  for  a chosen  point, 
a local  coordinate  system  is  defined  with  this  point  on  the  2-axis.  The  local  point  set 
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is  projected  onto  the  rcy-plane  and  the  surface  heights  provide  the  data  values.  This 
typically  gives  a configuration  very  close  to  the  hexagon  points  (3.1)  with  an  additional 
point  at  the  centre,  except  for  those  points  associated  with  the  icosahedron  vertices  where 
the  arrangement  is  a pentagon.  As  the  iscosahedral  mesh  is  refined  these  configurations 
become  less  regular. 

The  addition  of  the  central  point  to  the  hexagon  points  increases  the  order  of  the 
approximation.  When  the  surface  is  spherical,  the  symmetry  of  the  data  ensures  that  the 
computed  unit  normal  at  the  centre  point  for  polynomial,  MQ  or  TPS  is  exact  except  for 
rounding  error  ( e.g . for  MQ  the  error  is  ||n  — «mq||2  = 3 x 10~14).  However,  taking  MQ 
with  c = 10  and  a sphere  of  radius  9,  if  the  central  point  is  displaced  from  the  origin  to 
(0.01,0.01)  the  error  in  the  normal  is  3 x 10-3.  To  illustrate  convergence  for  an  irregular 
point  set,  the  hexagon  points  are  perturbed  by  the  addition  of  a factor  (i  — l)eh[l,  1]T 
for  points  i = 1, 2 . . . , 6 with  h the  radius  of  the  circumcircle  and  taking  e = 0.05.  For 
MQ  with  c = 10,  the  error  in  the  surface  normal  is  0(h 3)  whereas,  for  c = 0.4,  the  error 
is  larger  and  the  rate  of  convergence  varies  (see  Table  3). 


h 

||n-  nA/Q||2,  c=  10 

\\n-nMQ\\2,  c = 0.4 

1.0 

0.5 

0.1 

0.05 

3.15  x 10-b 
3.85  x 10~6 
3.06  x 10~8 
4.99  x 10“9 

6.14  x 10~3 
1.26  x 10~3 
6.35  x 10~6 
8.47  x 10-7 

Tab.  3.  Error  in  MQ  approximation  to  surface  normal  of  sphere,  irregular  point 
set. 

Accurate  curvature  values  are  essential  for  an  interface  which  is  driven  by  surface 
tension.  The  exact  value  of  k = -2/9  for  a sphere  of  radius  9,  together  with  the  computed 
values,  are  shown  in  Table  5.  The  polynomial  and  MQ  with  c = 10  are  close  to  the  exact 
value. 


Method 

K 

exact 

polynomial 
MQ  c = 0.1 
MQ  c = 10.0 

-0.222. . . 
-0.222912 
-1.638002 
-0.225387 

Tab.  4.  Curvature,  k = 2H , evaluated  at  the  central  point  of  a regular  hexagon. 

It  is  found  that,  for  the  icosahedral  mesh  with  N = 362,  the  local  point  sets  are 
sufficiently  regular  to  give  good  accuracy  for  surface  normals  and  curvature  using  MQ 
interpolants  when  c is  chosen  to  be  ‘large’  in  relation  to  the  point  spacing.  This  mesh  also 
gives  a corresponding  accuracy  for  the  discretised  integral  equation.  These  points  can 
thus  be  considered  ‘good’  for  the  MQ  approximation.  However,  if  the  mesh  is  further 
refined  or  the  surface  deforms  during  its  evolution,  then  the  approximation  becomes 
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‘less  good’  as  the  regularity  of  the  point  locations  is  lost.  Numerical  experiments  suggest 
second  order  convergence  with  point  separation  for  irregular  local  point  sets. 

6 Conclusions 

The  behaviour  of  MQ  and  TPS  interpolants  can  be  interpreted  by  reference  to  the 
corresponding  ‘least’  polynomial  interpolant,  with  the  MQ  connecting  the  polynomial 
C°°  surface  to  the  tensioned  surface  of  the  TPS  as  the  parameter  c decreases.  The  MQ 
interpolant  with  ‘large’  c (relative  to  the  point  separation)  exhibits  the  properties  of  the 
polynomial  case  and  is  similarly  affected  by  the  location  of  data  points.  Thus,  points 
on  a circle  in  the  plane  can  be  ‘good’  if  the  function  to  be  represented  is  harmonic, 
but  in  general  give  only  first  order  convergence  on  the  interior.  For  data  on  the  sphere, 
‘good’  points  for  polynomial  interpolation  are  also  good  for  the  MQ  with  ‘large’  c,  but 
other  near  equispaced  point  distributions  appear  to  give  similar  accuracy  with  MQ. 
The  tensioning  effect  of  smaller  values  of  c can  improve  the  results  if  the  underlying 
function  is  not  C°°.  When  applied  to  an  evolving  interface,  starting  from  an  initially 
spherical  shape  and  a refined  icosahedral  point  distribution,  it  is  found  that  local  MQ 
approximations  to  the  surface  derivatives  are  affected  by  the  point  locations.  This  can 
be  understood  by  reference  to  the  polynomial  interpolant  to  data  located  on  a circle  and 
causes  an  irregularity  in  the  convergence  as  N increases. 
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